neighborhoods

vizinhanças por contiguidade

Carregando shapefile de são paulo:

library(rgdal)

## Carregando pacotes exigidos: sp

## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.0, released 2018/12/14
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-5

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

shp_sp <- readOGR(dsn = "data/shapefiles/SP_Municipios", layer = "SP", verbose = F)
plot(shp_sp)

Critério queen de vizinhança:

library(spdep)

## Carregando pacotes exigidos: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Carregando pacotes exigidos: sf

## Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0

vizinhos_queen = poly2nb(pl = shp_sp, queen = TRUE,row.names = shp_sp@data$NM_MUN)

plot(shp_sp, border = "lightgray")
plot(vizinhos_queen, 
     coordinates(shp_sp), 
     add = TRUE, 
     col = "#33638DFF")

Matriz W Binária:

matrizW_queen <- nb2mat(neighbours = vizinhos_queen, style = "B", zero.policy = T)
colnames(matrizW_queen) <- shp_sp@data$NM_MUN

Gerando uma lista com várias matrizes W, cada uma com contiguidades de ordem 1 até 5:

vizinhos_queen_ordens <- nblag(neighbours = vizinhos_queen, maxlag = 5)

plot das vizinhanças com ordem 3:

plot(shp_sp, border = "lightgray")
plot(vizinhos_queen_ordens[[3]], 
     coordinates(shp_sp), 
     add = TRUE, 
     col = "#33638DFF")

Estabelecendo vizinhanças por contiguidade, critério rook (torre, não pega diagonais):

vizinhos_rook <- poly2nb(pl = shp_sp, queen = FALSE,row.names = shp_sp@data$NM_MUN)

plot(shp_sp, border = "lightgray")
plot(vizinhos_rook, 
     coordinates(shp_sp), 
     add = TRUE, 
     col = "#95D840FF")

Matriz W do critério rook:

matrizW_rook <- nb2mat(neighbours = vizinhos_rook, style = "B", zero.policy = TRUE)
colnames(matrizW_rook) <- shp_sp@data$NM_MUN

Vizinhanças por Distância Geográfica

Shapefile do estado da Bahia:

shp_ba <- readOGR(dsn = "data/shapefiles/shapefile_ba/", layer = "ba_state", encoding = "UTF-8", use_iconv = TRUE, verbose=F)

vizinhanças por distâncias geográficas (as distâncias serão calculadas em km):

library(spdep)
vizinhos_distancias <- dnearneigh(coordinates(shp_ba), d1 = 0, d2 = 90, longlat = TRUE)
summary(vizinhos_distancias)

## Neighbour list object:
## Number of regions: 417 
## Number of nonzero links: 11662 
## Percentage nonzero weights: 6.706577 
## Average number of links: 27.96643 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  5  6  6  2  4  4  7 10  6 16  6 10  7 10  4 12  8 17 18 20 15 11 12 11  7  5 
## 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 
##  2  3  4  9  5  3  8  7  2  8  5  5  5  5  2  7  4  8  6  3  8  5  3  1  4  8 
## 53 54 55 56 57 58 59 60 61 62 63 64 65 67 
##  3  4  3  3  3  5  8  2  3  4  5  3  1  1 
## 5 least connected regions:
## 113 393 408 409 415 with 1 link
## 1 most connected region:
## 240 with 67 links

Vizualizando as distâncias:

plot(shp_ba, border = "lightgray")
plot(vizinhos_distancias, coordinates(shp_ba), col = "#CC6A70FF", add =T)

Matriz W:

matrizW_distancias <- nb2mat(neighbours = vizinhos_distancias, style = "B")
colnames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
rownames(matrizW_distancias) <- shp_ba@data$MUNICIPIO

Vizinhanças Ponderadas por k-Nearest Neighbors

library(rgdal)

shp_sc = readOGR("data/shapefiles/shapefile_sc", layer="sc_state", verbose=F)

Obrigando cada municipio ter exatamente 3 vizinhos:

lista_knear <- knearneigh(coordinates(shp_sc), longlat = TRUE, k = 3)
vizinhos_knear <- knn2nb(knn = lista_knear)

plot(shp_sc, border = "lightgray")
plot(vizinhos_knear, coordinates(shp_sc), add = T, col = "#13306DFF")

Matriz Binária W:

matrizW_knear  <- nb2mat(neighbours = vizinhos_knear, style = "B")
colnames(matrizW_knear) <- shp_sc@data$NM_MUNICIP
rownames(matrizW_knear) <- shp_sc@data$NM_MUNICIP

Padronização em Linha da Matriz W

matrizW_queen_linha <- nb2mat(vizinhos_queen, style = "W", zero.policy = T)
colnames(matrizW_queen_linha) <- rownames(matrizW_queen_linha)

São paulo tem 23 vizinhos no método queen:

sum(matrizW_queen["São Paulo",])

## [1] 23

Quando padronizamos por linha obrigamos a soma da linha dar 1:

sum(matrizW_queen_linha["São Paulo",])

## [1] 1

Dupla Padronização da Matriz W

matrizW_queen_dupla_padr <- nb2mat(vizinhos_queen, style = "U", zero.policy = T)
colnames(matrizW_queen_dupla_padr) <- rownames(matrizW_queen_dupla_padr)

Quando padronizamos duplamente a soma por linha não da mais 1:

sum(matrizW_queen_dupla_padr["São Paulo",])

## [1] 0.006277293

Mas a soma da matriz inteira da um:

sum(matrizW_queen_dupla_padr)

## [1] 1

Padronização da Matriz W pela Estabilização da Variância

matrizW_queen_est_var <- nb2mat(vizinhos_queen, style = "S", zero.policy = T)
colnames(matrizW_queen_est_var) <- rownames(matrizW_queen_est_var)

Nesse caso a soma de matriz nos da a quantidade de municipios com vizinhos (ilha bela não tem vizinhos no método queen):

sum(matrizW_queen_est_var)

## [1] 644

Autocorrelação Espacial para Matriz W

Tipos:

• Auto correlação global
• AUto correlação local (granularidade no nível da observação)

Auto correlação global

Recuperando IDH de são Paulo e colocando no shapefile:

sp = read.csv("data/sp.csv")
shp_sp <- merge(x = shp_sp, y = sp, by.x = "CD_MUN", by.y = "codigo")

Para o cálculo da Estatística I de Moran, nosso algoritmo esperará como declaração um objeto de classe listw:

listw_queen <- mat2listw(matrizW_queen)

Estamos verificando a correlação de uma variável contra ela mesma defasada espacialmente (ou seja, comparando com os vizinhos) em relação a uma determinada observação:

• Hipótese Nula: Não existe auto correlação global, ou seja, os valores do IDH são independentes dos vizinhos
• Hipótese Alternativa: Há correlação. Valores do IDH são dependentes dos valores dos vizinhos

Teste da Hipótese: o valor cálculado de I (0.2328220224) é estatisticamente ao valor esperado de I (0.00155521)? Sim, eles são estatisticamente diferente e rejeitamos a hipótese nula.

Valor esperado de I :

-1/(644-1)

## [1] -0.00155521

Estatística I de moran:

idhs = shp_sp@data$idh
idhs_scale = scale(idhs)
I = (t(idhs_scale) %*% matrizW_queen %*%  idhs_scale /
    t(idhs_scale) %*% idhs_scale) *
  ( (length(idhs)-1)/sum(matrizW_queen) )
I

##          [,1]
## [1,] 0.232822

Como p-value = 6.695e-10 rejeitamos a hipótese nula e a estatsitica I de autocorrelação global é: 0.2328220224 - ou seja, existe autocorrelação global e ela é positiva.

moran.test(x = shp_sp@data$idh, listw = listw_queen, zero.policy = T)

## 
##  Moran I test under randomisation
## 
## data:  shp_sp@data$idh  
## weights: listw_queen  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 10.092, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.2328220224     -0.0015552100      0.0005393647

O Diagrama da Estatística I de Moran:

moran.plot(x = shp_sp@data$idh, listw = listw_queen, zero.policy = TRUE, 
           xlab = "IDH", ylab = "IDH Espacialmente Defasado")

Interpretação:

• As linhas que dividem os quadrantes são as médias de IDH em cada eixo
• Como estamos trabalhando com a matriz de contigunidade, ilha bela não tem viznho, por isso ele está destacado com o simbolo diferente no diagrama.
• São paulo está no quadrante High-High, ou seja, tem similaridade (transbordamento espacial) e tem valores alto da variável de interesse, IDH, com vizinho que também tem IDH alto
• àguas de são pedro está no quadrante High-Low e possui valores altos de IDH com vizinho com valores baixos
• Riversul está no quadrante Low-Low e tem similaridade (transbordamento espacial) e tem valores baixo da variável de interesse, IDH, com vizinhos que também tem IDH baixo
• Ibiuna está no quadrante Low-High e possui valores baixo de IDH com vizinhos com valores altos

Diagrama de Moran com ggplot:

idh_defasada = lag.listw(x=listw_queen, var=shp_sp@data$idh,zero.policy = T)

df = data.frame(cidade = shp_sp@data$NM_MUN, idh=shp_sp@data$idh,idh_defasada)
df["idh_defasada_na_unha"] = matrizW_queen %*% shp_sp@data$idh

plotly::ggplotly(
  df |>
    ggplot(aes(label=cidade)) +
    geom_point(aes(x=idh,y=idh_defasada), alpha=0.5,size=3) +
    geom_smooth(aes(x=idh,y=idh_defasada), method="lm", se=F) +
    geom_hline(yintercept = mean(df$idh_defasada),lty=2) +
    geom_vline(xintercept = mean(df$idh),lty=2) +
    annotate("text", x=0.66, y=18.5, label="Low-High") +
    annotate("text", x=0.84, y=18.5, label="High-High") + 
    annotate("text", x=0.66, y=-1, label="Low-Low") + 
    annotate("text", x=0.84, y=-1, label="High-Low") + 
    theme_bw()
)

## `geom_smooth()` using formula 'y ~ x'

Auto correlação local

1. devemos padronizar a matriz W em linha:

matrizW_queen_linha <- nb2mat(vizinhos_queen, style = "W",zero.policy = T)
colnames(matrizW_queen_linha) = rownames(matrizW_queen_linha)
listw_queen <- mat2listw(matrizW_queen_linha)

Estatística Moran Local com o uso da função localmoran:

moran_local <- localmoran(x = shp_sp@data$idh, listw = listw_queen, zero.policy = T)

• Cluster Espaciais: O pontos com Moran Local estatisticamente significantes e que estejam em High-High ou Low-Low
• Outluers Espacial: O pontos com Moran Local estatisticamente significantes e que estejam em High-Low ou Low-High

idh_zscore = scale(shp_sp@data$idh)
moran_local_na_unha = rowSums(sweep(x=matrizW_queen_linha,MARGIN=2,STATS = idh_zscore, FUN = "*"))*idh_zscore

Juntando os resultados da Estatística Moran Local no dataset do objeto shp_sp:

moran_local_mapa <- cbind(shp_sp, moran_local)

library(gtools)
library(tmap)
quantile(moran_local_mapa@data$Ii, type = 5, probs = c(0,0.2,0.4,0.6,0.8))

##           0%          20%          40%          60%          80% 
## -1.424372672 -0.123143140  0.002717804  0.121068801  0.458396014

moran_local_mapa@data <- moran_local_mapa@data %>% 
  mutate(faixa_quantis = factor(quantcut(x = Ii, q = 5))) 

tm_shape(shp = moran_local_mapa) +
  tm_fill(col = "faixa_quantis", palette = "-magma") +
  tm_borders()

Considerando os quadrantes:

moran_local_mapa <- cbind(moran_local_mapa, 
                          attr(x = moran_local, which = "quadr")[1])

tm_shape(shp = moran_local_mapa) +
  tm_fill(col = "mean", palette = "-viridis") +
  tm_borders(col = "gray")

Criando um vetor que contenha o centro das observações da variável idh ao redor de sua média:

quadrantes <- vector(mode = "numeric", length = nrow(moran_local))
centro_idh <- shp_sp@data$idh - mean(shp_sp@data$idh)
centro_moran_local <- moran_local[,1] - mean(moran_local[,1])

Vamos manter na base somente os I locais significantes. Enquadrando nossas observações em seus respectivos quadrantes:

sig <- 0.05
quadrantes[centro_idh > 0 & centro_moran_local > 0] <- "HH"
quadrantes[centro_idh > 0 & centro_moran_local < 0] <- "HL"
quadrantes[centro_idh < 0 & centro_moran_local > 0] <- "LH"
quadrantes[centro_idh < 0 & centro_moran_local < 0] <- "LL"

Ajustando a presença da observação em razão de sua significância estatística:

quadrantes[moran_local[,5] > sig] <- "Estatisticamente_não_significante"
moran_local_mapa@data["quadrantes"] <- factor(quadrantes)

Plotando os quadrantes de forma espacial (versão ‘default’):

tm_shape(shp = moran_local_mapa) +
  tm_polygons(col = "quadrantes", 
              pal = c(HH = "darkred",
                      HL = "red", 
                      LH = "lightblue", 
                      LL = "darkblue",
                      Estatisticamente_não_significante = "white")) +
  tm_borders()

## Warning: One tm layer group has duplicated layer types, which are omitted. To
## draw multiple layers of the same type, use multiple layer groups (i.e. specify
## tm_shape prior to each of them).

Autocorrelação Local - A Estatística G de Getis e Ord

Recarregando dados da Bahia:

library(rgdal)
library(tidyverse)
library(spdep)

shp_ba <- readOGR(dsn = "data/shapefiles/shapefile_ba/", layer = "ba_state", 
                  encoding = "UTF-8", use_iconv = TRUE, verbose=F)
vizinhos_distancias <- dnearneigh(coordinates(shp_ba), d1 = 0, d2 = 90, longlat = TRUE)
matrizW_distancias <- nb2mat(neighbours = vizinhos_distancias, style = "B")

colnames(matrizW_distancias) <- shp_ba@data$MUNICIPIO
rownames(matrizW_distancias) <- shp_ba@data$MUNICIPIO

idh_bahia = read.csv("data/idh_bahia.csv")

Adicionando dados de IDH da bahia ao shapefile:

idh_bahia$Codigo = as.character(idh_bahia$Codigo)
shp_ba@data = shp_ba@data |> left_join(idh_bahia, by = "Codigo")

Transformação do objeto matrizW_distancias em um objeto de classe listw:

listw_dist <- mat2listw(x = matrizW_distancias)

Calculando a Estatística G de Getis e Ord:

g_local <- localG(x = shp_ba@data$idh, listw = listw_dist)

Juntando as informações do objeto g_local ao nosso shapefile:

mapa_G <- cbind(shp_ba, as.matrix(g_local)) 
mapa_G@data = mapa_G@data |> rename(estistica_g = 4)

Plotando a Estatística G de forma espacial:

library(tmap)
tm_shape(mapa_G) + 
  tm_fill("estistica_g", palette = "-RdBu") + 
  tm_borders()

Estradificação por 8 e melhorado para daltônicos:

library(gtools)
mapa_G@data = mapa_G@data |> mutate(faixa_quantis = factor(quantcut(x = estistica_g, q = 8))) 

tm_shape(mapa_G) + 
  tm_fill("faixa_quantis", palette = "plasma") + 
  tm_borders()